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We present a Bayesian hierarchical modeling approach to pa- 
leoclimate reconstruction using borehole temperature profiles. The 
approach relies on modeling heat conduction in solids via the heat 
equation with step function, surface boundary conditions. Our anal- 
ysis includes model error and assumes that the boundary conditions 
are random processes. The formulation also enables separation of 
Q, | measurement error and model error. We apply the analysis to data 

^* ■ from nine borehole temperature records from the San Rafael region 

in Utah. We produce ground surface temperature histories with un- 
certainty estimates for the past 400 years. We pay special attention 
to use of prior parameter models that illustrate borrowing strength 
in a combined analysis for all nine boreholes. In addition, we review 
selected sensitivity analyses. 



^~) 1. Introduction. Reconstruction of past climate plays an important role 

in climate change analysis. Comparisons between climate behavior before 
■r^J- ■ and after human influences are a relevant component of claims of attribution 

r~^ . of climate change to our activities. The term paleoclimate is used in reference 

to data analysis and modeling of climate for times before the modern era of 
data collection. The time periods of interest range from hundreds to millions 
of years before the present. Of course, as our interest moves toward the past, 
the availability of reliable and spatially and temporally plentiful observations 
of weather and climate diminishes. In response, scientists have developed the 
use of proxy indicators of climate [Jansen et al. (2007)]. 

A proxy is a quantity taking values that respond to climate behavior. For 
example, annual tree ring thicknesses respond to weather variables that con- 
trol growth, that is, temperature and precipitation. If one develops useful 
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models for proxies as rough functions of climate, then the models can be in- 
verted to estimate climate behavior based on observed proxy variables. The 
statistical notions of regression and inverse regression analyses are immedi- 
ately evident. Though exceptions exist and the trend is positive [e.g., Haslett 
et al. (2006); Li, Nychka and Ammann (2007, 2010)], there has been insuf- 
ficient participation by statisticians in paleoclimate reconstruction. This is 
surprising in view of the richness of the statistical challenges: proxies are 
themselves observed with error; the forward models for proxies as func- 
tions of climate are partially known at best and subject to model errors; 
inverse analyses are not trivial statistically; spatial and temporal coverage 
and mismatches between proxy data sets and desired climate inferences 
are among the issues. Furthermore, there is substantial interest in paleo- 
climate among policy makers and the general public [Wegman, Scott and 
Said (2006); Smith, Berliner and Guttorp (2010)]. 

In this article we focus on the critical problem of surface temperature re- 
construction [Jansen et al. (2007); North et al. (2006)]. We analyze borehole 
temperature data sets and their use in reconstructing surface temperature 
time series [e.g., Beltrami and Mareschal (1995); Pollack, Huang and Shen 
(1998)]. A borehole is a narrow shaft drilled into the ground (or ice), typi- 
cally vertically, in search of subterranean resources (gas, oil, water, minerals, 
etc.). Boreholes are also used to monitor environmental processes (e.g., per- 
colation of contaminants) or as pilots to access suitability for more intense 
drilling or construction projects. Borehole data that are used for tempera- 
ture reconstruction are typically obtained as byproducts of such projects. 
Therefore, borehole data are observations of opportunity rather than having 
been designed with climate reconstruction in mind. We note that borehole 
data are temperature measurements, and, hence, perhaps not as indirect 
a measure of surface temperature as other proxies. However, the problems 
of developing and inverting a model for borehole data as a function of surface 
temperatures are challenging. 

The underlying theory for using borehole data to infer surface tempera- 
ture is the physics of heat conduction. In principle, the transfer of heat is 
governed by the heat equation. This is a partial differential equation describ- 
ing the temporal evolution of the temperature field over some domain. The 
idea is that the surface temperatures over time serve as boundary condi- 
tions for the evolution of temperature below the surface. Then information 
regarding subsurface temperatures can be inverted to estimate the boundary 
conditions. 

There are important issues and uncertainties that arise in applying this 
strategy. First, subsurface temperatures respond to ground-surface tempera- 
tures as opposed to near surface air temperatures. Though the latter two are 
related, they are not identical, perhaps due to snow cover and other factors. 
Next, as heat conducts into deepening levels of the subsurface, it spreads or 



TEMPERATURE RECONSTRUCTION 3 

smears out, leading to losses in information regarding the boundary as time 
increases. Though this problem is well known, we will seek explicit charac- 
terizations of this loss of information as reflected in uncertainty measures 
associated with our results. Another issue is that there are factors affecting 
heat conduction that are difficult to quantify. For example, conduction rates 
depend on characteristics of the media (i.e., rock types) through which the 
heat flows. Further, percolation of water through the media also impacts 
heat flow. For such reasons, we incorporate the heat equation with error 
and unknown parameters in our modeling. 

We present Bayesian modeling and analysis for data from 9 boreholes in 
the San Rafael region in Utah. To combine information from these boreholes, 
we assume model parameters are site-specific, but sampled from common 
distributions. Our modeling incorporates both the observations and physics 
into an analysis that is sensitive to the uncertainties in both information 
sources. The use of such physical-statistical analyses in geophysical problems 
is increasing; for examples, see Berliner (2003), Berliner et al. (2008) and 
Wikle et al. (2001). See Hopcroft, Gallagher and Pain (2007) for a related 
Bayesian analysis of borehole data. 

1.1. Review: Borehole data analysis. The conventional approach is to 
frame analyses in terms of reduced temperatures defined as follows. For 
a given borehole, consider N depths zi,...,zn, where increasing values of z 
correspond to increasing depths. Let T be the N x 1 vector of true temper- 
atures at these depths. The corresponding vector of reduced temperatures 
is given by 

(1) T r = T - T Q 1 - q K, 

where To is the surface temperature intercept, 1 is an N x 1 vector whose 
elements are all equal to one, qo represents background heat flow, and R is 
anJVxl vector of thermal resistances at each of the depths. This modeling 
step is intended to account for the fact that both heating from the earth's 
core and rates of heat conduction vary with depth, thereby justifying use of 
the simple heat equation model described in Section 3.1. 

The thermal resistances account for differences in heat conduction and are 
assumed known throughout the analysis (see Section 2). To deal with the 
unknowns To and qo, it is customary to replace the true temperatures T by 
the observed temperatures, Y, in (1). Then, To and qo are estimated via least 
squares by regressing Y onto R. In that step a subset of the data is used, 
corresponding to those depths where the climate change signal is assumed 
to be negligible [for our data this means below 150 m or 200 m, depending 
on the region; Harris and Chapman (1995)]. The resulting estimates 

(2) f r (z i )=Y(z i )-(f +q R(z l )), i = l,...,N, 
are then treated as the true reduced temperatures. 
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Having made the above adjustments, the heat equation is assumed to 
apply. Let T/j be a K x 1 surface temperature history vector. Here, the 
surface temperatures are assumed to be constants over K time intervals used 
in the analysis. The heat equation can be solved (see Section 3), leading to 
the linear relationship 

(3) t r = AT h , 

where A is an N x K matrix developed from the solution to the heat equation 
[see (10)]. The objective then is to solve the inverse problem, that is, obtain 
an estimate of TV In most examples, A is ill-conditioned and the inversion 
of A' A is unstable. Therefore, traditional regression methods which involve 
taking the inverse of A' A lead to unstable estimates of T/j. A common 
approach is to use a singular value decomposition (SVD) of the A matrix 
and retain only a few of the singular vectors [Vasseur et al. (1983); Beltrami 
and Mareschal (1991); Mareschal and Beltrami (1992); Harris and Chapman 
(1995)]. In this paper we take a hierarchical Bayesian regression approach 
that does not involve taking the inverse of A' A. Hence, we avoid having to 
pick the number of singular vectors (principal components) to retain. 

Other approaches can be found in the geophysical literature. For example, 
functional space inversion is a popular method [Shen and Beck (1991, 1992); 
Harris and Chapman (1998)]. A comparative study of some inverse methods 
in this setting can be found in Shen et al. (1992). 

1.2. Outline. The paper is organized as follows. In Section 2 we describe 
the borehole data used in the analysis. In Section 3 after a brief introduction 
to the physical model the analysis is based on, we develop a Bayesian hierar- 
chical model. To best convey the ideas, we first present a single-site borehole 
model in Section 3.2. We extend it to include data from multiple boreholes 
in Section 3.3. The results are presented in Section 4. In Section 4.2 we com- 
pare the results to those of single-site models and in Section 4.3 we present 
a number of sensitivity analysis. We end with a discussion in Section 5. 

2. Data. We consider borehole data from the Colorado Plateau in Utah. 
The data consist of nine measured temperature-depth profiles (shown in Fig- 
ure 1) belonging to two regions, the San Rafael Desert and the San Rafael 
Swell. The geography of these regions is characterized by layered sedimen- 
tary rocks that each have different thermal conductivities. Measurements 
of these different conductivities are available [Bodell and Chapman (1982)] 
and have been adjusted to the specific formations in the boreholes used in 
this analysis so that the estimated thermal conductivity for each forma- 
tion may be different between regions but not within regions [see Harris 
and Chapman (1995) for details about these adjustments]. The adjusted 
thermal conductivities (k) for each sedimentary formation and abbreviated 
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San Rafael Desert 



San Rafael Swell 



14.51 15.73 16.00 16.17 



O- 11.3 12.08 12.54 13.54 



Q 



. SRD-3 SRD-4 
SRD-2 



SRD-1 

Relative temperature [ 9 C] 




Relative temperature [ e C] 



Fig. 1. Measured temperature- depth profiles from boreholes in the San Rafael Desert 
(left) and the San Rafael Swell (right) . The temperatures are shifted so they do not overlap, 
one tick on the x axis corresponds to 1° C. The value of the shallowest measurement is 
shown above each profile. The horizontal line segments show the formation boundaries, the 
names of the formations are given in Table 1. 



formation names are shown in Table 1 along with the formation boundaries 
within each borehole. For definitions of the abbreviated formation names 
see Bodell and Chapman (1982). The sedimentary formation boundaries are 
also shown as small horizontal line segments in Figure 1. For more back- 
ground on the data and temperature reconstructions based on them see 
Harris and Chapman (1995, 1998). 

Thermal resistance (R) is a function of depth and the thermal conductiv- 
ity at that depth, and is calculated as 



(4) 



R(zi 



E 
i=i 



zi - Zl-i 



k(zi] 



1. 



.,N, 



where k{zi) is the conductivity for the depth interval from z\^\ to Z[. Recall 
that zq = denotes the surface and z is increasing in depth. If the conduc- 
tivity is constant (k) for the entire borehole, (4) reduces to R{zi) = Zi/k, 
that is, a constant temperature gradient. 

3. Bayesian hierarchical model. 



3.1. Physics based modeling. Following convention, we assume that for 
reduced temperatures, boreholes are well approximated as homogeneous, 
heat source free (except at the surface), one-dimensional, semi-infinite solids 
(i.e., a single boundary is at the surface). It follows that the reduced tem- 
peratures, T r {z,t) at depth z and time t, can be reasonably modeled by the 
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Table 1 

Depth to formation boundaries within each borehole, names of the 

formation types down to each boundary and the corresponding thermal 

conductivities (k). The table also shows for each borehole the estimated 

surface temperature intercept, To, and the year the borehole was 



Borehole 


Depth [m] 


Formation 


fe [W/mK] 


San Rafael Desert 








SRD-1 


60 


Jca 


2.91 


fa = 13.72 


225 


Jna 


4.09 


year: 1979 


260 


JTrk 


3.96 




395 


Trwi 


3.86 


SRD-2 


25 


Jca 


2.91 


To = 15.12 


215 


Jna 


4.09 


year: 1976 


275 


JTrk 


3.96 




365 


Trwi 


3.86 


SRD-3 


140 


Jna 


4.09 


To = 15.38 


200 


JTrk 


3.96 


year: 1979 


320 


Trwi 


3.86 


SRD-4 


145 


Jna 


4.09 


To = 15.51 


210 


JTrk 


3.96 


year: 1979 


320 


Trwi 


3.86 


SRD-7 


185 


Jna 


4.09 


To = 13.16 


260 


JTrk 


3.96 


year: 1980 


375 


Trwi 


3.86 


San Rafael Swell 








SRS-3 


250 


Pco 


5.01 


To = 10.76 


390 


Pec 


4.35 


year: 1979 


400 


Mr 


4.82 


SRS-4 


135 


Jca 


2.91 


To = 11.82 


375 


Jna 


4.18 


year: 1979 


410 


JTrk 


3.86 




510 


Trwi 


4.17 


SRS-5 


55 


Jca 


2.91 


To = 11.82 


350 


Jna 


4.18 


year: 1979 


400 


JTrk 


3.86 




480 


Trwi 


4.17 


WSR-1 


50 


.Is 


4.10 


% = 12.87 


105 


Jcu 


3.96 


year: 1980 


245 


Je 


3.43 




320 


Jca 


2.91 




455 


Jna 


4.18 




515 


JTrk 


3.86 




575 


Trwi 


4.17 
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heat equation, 

8 2 T r (z,t) 1 8T r (z,t) 

[ j a? 2 K at ' 

where k is the thermal diffusivity of rock. As is customary in borehole anal- 
ysis, we fix k to be 10 -6 m 2 /s [see, e.g., Harris and Chapman (1995)]. The 
boundary condition, T r (0, t), is the primary target of our inference. Assum- 
ing that the initial reduced temperatures, T r (z,t\), are zero for all depths z, 
the solution to (5) is 

2 r° „ / z 2 \ _„2 



(6) ^'^^l^H^-^r^ 

[e.g., Carlslaw and Jaeger (1959)]. We assume that the boundary function 
is a step function, 

Ti, ift 1 <t<t 2 , 

, N T 2 , if t 2 <t<t 3 , 

(7) T r (0,t) = < . 

,Tk, if tK <t< tR+l- 

Then the solution (6) reduces to 

r r (^,t) = r 1 erfc( r I > )+(T 2 -r 1 )erfcf ^ , ) 

(8) 

+ ... + ( r,-r,_ 1 ,e*(- 7 =i=), 

where erfc(-) is the complementary error function 
(9) erfc(x) = -= / e"^ da. 



x 

Here, the time points ti, ■ ■ ■ ,tx+i are calendar years, t\ being the earliest 
year considered and tx+i is the year the borehole data were collected. For 
example, borehole SRD-7 was logged in t±2 = 1980 and we selected the fol- 
lowing years for the time points ti,...,tu: 1600, 1650, 1700, 1750, 1800, 
1850, 1875, 1900, 1925, 1950 and 1965. This range is in concert with other 
analyses. We also performed some analyses using a slightly more refined 
temporal grid. These resulted in little change in the general form of the 
posterior results. We note that for highly refined grids, the structure and 
dimension of A becomes an issue. 

Collecting terms in (8), we can write the solution at t = tjc+i in vector 
form: 

(10) T r = AT h , 
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where T r = (T r (zi,t K+1 ), . . . ,T r (z N ,t K+1 ))' , T h = (T x , . . .,T K )', and A is 
an N x K matrix with (z,j)th entry 



erfc 



I — ; I ~~ er fc I 



\yjAK{t K +\ ~tj)J \y/AK(t K +l -tj+i] 

for i = 1, . . . , N and j = 1, . . . , K — 1 and (i, K)th entry 

erfcf Zi ); i = l,...,iV. 

\V 4 «(*Jf+i -*«■)/ 

3.2. Single-site model. It is useful to view the modeling in three basic 
stages, a data model, process model and a parameter model. 

(i) Data model. Let Y be a vector of observed temperatures at depths 
Z\ , . . . , zn and let T denote the iV-dimensional vector of corresponding true 
temperatures. We assume that the observations are noisy, unbiased mea- 
surements of the true temperatures. Specifically, we assume that 

(11) Y = T + e, 

where e is an iV-dimensional vector of normally distributed, independent 
errors, all with mean zero and common variance a\ (see Section 4.3.4 for 
discussion of the independence assumption). 

Recalling the definition of reduced temperatures in (1), the true temper- 
atures can be written as 

(12) T = T r + T l 7V + (?oR, 

where T r is the vector of true reduced temperatures and other quantities 
are defined after (1). 

Combining (11) and (12), the assumed data model is the conditional dis- 
tribution 

(13) Y\T r , q ,a^ ~ N(T r + T 1 N + q R,a^I N ), 

where the vertical bar | is read "given," ~ is read "is distributed," and In 
is the N x N identity matrix. 

Note that To and T r in (13) are not identifiable in that, for any constant c, 
the shifted parameters T r — > T r + cljv, To — > Tq — c yield identical data 
models. To circumvent this issue, we assume that To is known. The assumed 
values of To for the boreholes analyzed here are given in Table 1. We report 
on the sensitivity of results to the choice of To in Section 4.3.1. 

(ii) Process model. We let T^ denote the iC-vector containing the surface 
temperature history. We incorporate the heat equation in defining a stochas- 
tic process model 

(14) T r |T h ,£~JV(,4T h ,£) 
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[recall (10)]. We also assume a Gaussian prior for the histories: 

(15) T h ~N(n,T). 

(iii) Parameter model. We assume that the covariance matrix of the pro- 
cess model errors [see (14)] is diagonal with common variance, namely, 
£ = o~ 2 In- That is, after accounting for the dependence on the temperature 
history, the heat equation offers reliable explanation of reduced tempera- 
tures requiring only some local-in-depth errors. In Section 4.3.4 we consider 
sensitivity of results with respect to the independence assumption. 

Additional specifications of parameter priors is delayed until we discuss 
modeling for multiple sites. 

3.3. Multiple-site model: Spatially distributed parameters. We extend the 
model to combine data from multiple boreholes by allowing site-specific 
processes and parameters. A list of all the model parameters used in this 
model is provided in Appendix A. 

(i) Data model. We assume that measurements from different sites are 
conditionally independent with data models that depend only on site-specific 
processes and parameters. We also assume that reduced temperature vectors 
from different sites are conditionally independent with priors that depend 
on site-specific histories and parameters. Formally, the data model is the 
product of densities corresponding to the models 

(16) Yj | T rj , q oj , ay. ~ N N . (T rj + T 0j l Nj + qojRj , a\. I Nj ) 

for the 9 sites labeled j = 1, . . . , 9 with observation vectors of length Nj. 

(ii) Process model. Similarly, the process model for the reduced tempera- 
ture vectors is the product of densities for 

(17) T rj | T hj , a) ~ N Nj (AjT hj , ap Nj ) . 

Since all 9 sites are in the Colorado Plateau, we expect them to have been 
influenced by common large-scale climate effects. However, the sites are lo- 
cated in two subregions: Sites 1-5 are in the San Rafael Desert (D), Sites 6-9 
are in the San Rafael Swell (S). To account for common influences within 
subregions, we assume that all 9 histories are conditionally independent, 
with parameters depending on subregions: 

(18) T hj \ii D , 7 2 D ~N K (n D ,ry 2 D I K ), j = l,...,5, 
and 

(19) T hj \n s ,>y 2 s ~N K (vL S ,>y 2 s I K ), j = 6,...,9. 

We remark that this is a very elementary spatial model. More complex 
spatial modeling of parameters is feasible and recommended, depending on 
prior information and data richness. 
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(iii) Parameter model. 

Model for heat flow parameters: To account for region-wide influences, we 
assume that the heat flow parameters q = (qoi, • • • , qo9)' are sampled from 
priors with dependence structures. These priors are similar to exchangeable 
models [e.g., Section 4.6.2 in Berger (1985)]. 

The heat flow parameters are assumed to be conditionally independent 
with Gaussian priors where both the mean and the variance depend on which 
region the borehole is in, 

(20) qOj\VD,TD~N(v D ,Ti)), j = l,...,5, 
and 

(21) qoj\t>S,Ts~N(us,T%), i = 6,...,9. 

The means uq and us are assumed to be a priori independent with Gaussian 
priors having common mean v: 

(22) v D \v~ N(v,ri 2 ) and v s \v ~ N(is,r] 2 ). 
Next, we assume that 

(23) v~N(yo,r$). 

We can combine the distributions in (22) and (23) and integrate out v, 
leading to the following prior for vo and vs- 

™ (sH((:M"M"° /u 

Model for means of histories: We assume the means of the tempera- 
ture histories fi D and n s in (18) and (19) have common mean /x; /x D |/x~ 
N{h,o~ 2 d Ik) and fXg\fX ~ N([/,,o- s Ik)- We in turn assume that fi ^ N(fj, , 
UqIk)- As in the development of (24), integrating out /x yields the following 
joint prior for fi D and n s : 

Note that the covariance structures in (18), (19) and (25) only include some 
spatial dependence, but no temporal structure. Some spatial-temporal de- 
pendence among historical temperature values are displayed in their poste- 
rior distribution. 

Priors for the variances: The measurement error and process model er- 
ror variances appearing in (16)-(21) are all assigned independent, inverse 
gamma priors: 

(26) a Y , ~ IG(a Y ,b Y ) and a) ~ IG(a, b) for j = 1, . . . , 9, 



'j 



(27) Tf)~IG(a T ,b T ) and t§~ IG(a T ,b T ), 

l-p Oy). 



(28) j D ~ IG(a 7 ,6 7 ) and 7 S ~/G(a 7 ,6^ 
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3.4. Selection of parameters of prior distributions. We describe the se- 
lections of parameters of priors or hyper parameters introduced above: 

Measurement error variances (ay, by)- As suggested in Harris and Chap- 
man (1995), "The precision and accuracy of the measurements are estimated 
to be better than 0.01 K and 0.1 K, respectively." We view this as suggest- 
ing that a reasonable prior mean for the variances of the measurement er- 
rors, a Y ., j = 1, . . . ,9, is 0.11 2 . A conservative choice for the prior variances 
is 1.0. The values ay = 2.000146 and by = 0.012102 yield an inverse gamma 
distribution matching these properties. For additional intuition we remark 
that the 0.025 and 0.975 quantiles of this prior are equal to 0.002172 and 
0.049955, respectively. Further, the corresponding 0.025 and 0.975 quantiles 
for the ay. are 0.0466 and 0.2235, respectively. 

Model error variances (a, b). Though we know comparatively little about 
the variances a 2 of the model errors, we can develop some plausible expecta- 
tions. For example, if the standard deviations of the model errors are 0.5, we 
expect the model to be within 1.5°C from the truth 99.7% of the time. Hence, 
we specified the prior mean of each a 2 to be 0.50 2 and a very large prior vari- 
ance of 100. These selections correspond to a = 2.000625 and b = 0.250156. 
The corresponding 0.025 and 0.975 quantiles for the Gj are 0.212 and 1.016, 
respectively. 

Heat flow parameters (vQ,ii 2 D ,ri 2 s ,rj 2 l ,a T ,b T ). The background heat flow qo 
has been shown in other studies to range from about 30 mW/m 2 (milli- 
watt per meter-squared) to about 100 mW/m 2 , with the majority of values 
ranging between 50 and 70 mW/m 2 [Bodell and Chapman (1982); Beltrami 
and Mareschal (1995); Dorofeeva, Shen and Shapova (2002), e.g.]. Focus- 
ing on (24), we selected the prior mean vq = 60 mW/m 2 and set the stan- 
dard deviation 770 = 20 mW/m 2 . The standard deviations ?]d and rjs rep- 
resent variability due to subregion. We set r)n =775 = 10. Note that these 
selections imply that the prior standard deviations of vjj and v$ are equal 
to (20 2 + 10 2 ) - 50 m 22.36 mW/m 2 and the correlation between up and us 
is 0.80. We discuss sensitivities of results to these selections in Section 4.3.2. 

Recalling (20) and (21), t^, and rj quantify variability of the t/oj about 
their regionally defined prior means. We set the prior means for these vari- 
ances to 0.1 2 with a corresponding large prior variance of 1. It follows that we 
select a T = 2.000100 and b T = 0.010001. Note that the units here are W/m 2 , 
so this corresponds to r^, and r| having prior mean of 100 2 (mW/m 2 ) 2 
and the 0.025 and 0.975 quantiles for Try and ts are 42.4 mW/m 2 and 
203.2 mW/m 2 , respectively. 

Histories (/x , cr 2 D , a s , a 2 ,, a 7 , 6 7 ). In the model parameterizations used here 
both reduced temperatures and temperature histories represent departures 
from the baseline surface temperature To [see (12)]. Hence, a reasonable 
prior mean for T^ is fj, = 0. Further, these departures occur over time 
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intervals of lengths between 10 and 50 years. Focusing on (25), we selected 
ctq = 0.1 and a 2 D = cr| = 0.2. These selections imply that the prior standard 
deviations of the coordinated nr>k and nski k = l,...,K, are equal to (0.1 + 
0.2) ' 50 ~ 0.5477 and the correlation between /Ud& and nsk is 0.333. We 
discuss sensitivities of results to these selections in Section 4.3.3. 

Finally, recalling (18) and (19), 7^ and 7^ quantify variability of the 
elements of the history vectors about their regionally defined prior means. 
We set the prior means for these variances to 0.8 and prior variances equal to 
1. It follows that we select a 7 = 2.064 and 6 7 = 0.8512 and the corresponding 
0.025 and 0.975 quantiles for 7^ and 7^ are 0.387 and 1.801, respectively. 

4. Results. 

4.1. Multiple-site model. The multiple-site model described in Section 3.3 
has 818 unknown parameters, including 664 elements of the reduced temper- 
ature vectors T r j . We implemented a Gibbs sampler in R to obtain samples 
from the posterior distribution. The full conditional distributions used are 
given in Appendix B. We obtained 30,000 samples and discarded the first 
2000 as burn-in, leaving 28,000 samples to use for inference. Trace-plots for 
parameters and a random selection of elements of T r j showed no indication 
of convergence problems. We estimated marginal posterior densities using 
(Gaussian) kernel density estimation. 

The ground surface temperature (GST) histories, Tjy, are of primary 
interest. Estimated posterior means and credible sets for the GST histories 
are shown in Figure 2. In Figure 3 we show 5 samples of T^j for each 
borehole. These 5 samples are taken 6000 MCMC iterations apart (after 
burn-in) and 100 iterations apart between boreholes. Note that the generated 
realizations are not overly smooth, but that the spreads of the realizations 
are decreasing as time approaches the present. Estimated posterior means 
and credible sets for the mean GST histories for the San Rafael (SR) Desert 
and SR Swell regions, fx D , fj, s , and prior and posterior densities of the 
parameters jd and 7s, are shown in Figure 4. One striking feature of the 
GST histories is that posterior uncertainties are substantial. We see patterns 
of warming in the last century (and cooling for site WSR-1), but there are 
large posterior uncertainties associated with the estimated trends at each 
borehole. On the other hand, the posterior uncertainty is lower for more 
recent times (especially at the last time point). 

We note that the temperature trends (posterior means) are similar within 
the SR Desert sites, but quite variable within the SR Swell region, most no- 
tably at sites SRS-5 and WSR-1. Furthermore, the posterior credible inter- 
vals are wider for boreholes in the SR Swell. This difference is also apparent 
in the density estimates of the standard deviation of the GST histories: 75 
is larger than ju (see Figure 4, right). The mean GST histories fi D and n s 
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Fig. 2. Estimated posterior means and credible sets for the ground surface temperature 
(GST) histories, Thj ■ The white squares show the posterior means of Thj and the vertical 
bars show symmetric 50% (thicker and black) and 90% (thinner and grey) posterior credible 
intervals. 

show slightly different temperature trends (Figure 4, left), but the posterior 
uncertainty is somewhat large and increasing as we go further back in time. 
Estimated marginal posterior densities of the site-wise measurement and 
model error standard deviations, ayj and <7j, are shown in Figure 5. The 
locations of these densities are quite different from the prior means. The ayj 
are of the order 0.03-0.05°C compared to their prior mean 0.11°C. The <%,- 
are of the order 0.05-0.15°C compared to the prior mean 0.5°C. An inter- 
esting pattern emerges in Figure 5. Both ayj and Gj have higher posterior 
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Fig. 3. Ensembles of ground surface temperature (GST) histories generated from the 
posterior distribution. For each borehole, the thick black lines show the posterior means of 
Thj and the five grey lines are five different samples of T;^ . 

uncertainty for boreholes in the SR Desert than the SR Swell. Note that the 
boreholes in the SR Swell have more measurements than those in the SR 
Desert (see Figure 1). The SRD-2 borehole has by far the fewest measure- 
ments and, as indicated in Figure 5, densities for ayj and (Tj for that site 
are wider and are closer to the prior than the other densities. 

Estimated posterior densities of the background heat flow qoj are shown 
in Figure 6. Posterior means and 90% credible intervals for qoj and the 
means vd and v$ are shown in Table 2. Estimated posterior densities of 
the means and standard deviations of the heat flow, z/£>, us, td and ts, 
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are shown in Figure 7. It is clear from Figure 6 that the background heat 
flow is lower for boreholes in the SR Desert than in the SR Swell (except 
for sites SRD-1 and SRS-3). The 9 posterior densities show varying degrees 
of posterior uncertainty, in large part in response to the amount of data in 
each borehole. The high values of the standard deviations T£> and ts (see 
Figure 7) indicate the wide range of the heat flow q^j within the regions. 

There is considerable posterior uncertainty regarding the mean heat flows 
V£> and vg (see Figure 7), especially when compared to the relatively pre- 
cise posterior densities for the 9 individual heat flows qoj. This is what we 
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Fig. 6. Estimated posterior densities of the heat flow qoj for boreholes from the San 
Rafael Desert (solid lines) and San Rafael Swell (long dashes). 

expect since there is substantial variation of the locations of these precise, 
individual densities. Another point is that the posterior means of ujj and v$ 
are surprisingly alike, and do not seem to correspond to the difference in 
heat flows for the two regions that is apparent in Figure 6. For example, the 
posterior mean of ujj is higher than all the posterior means of qoj in the SR 
Desert (see Table 2). We explain this behavior in Section 4.3.2. 

4.2. Comparison to single-site models. In our primary analysis we com- 
bined data from all boreholes in one hierarchical multiple-site model. The 
temperature histories for boreholes in the same region share the same mean 
and variance. Similarly, the heat flow parameters for boreholes within the 
same region share the same mean and variance. One rationale for doing this 
is that it enables us to learn about region-wide mean temperature histories 
and region-wide mean heat flow. Another rationale is that hierarchically 
linking boreholes within regions allows for sharing of information between 
boreholes through the shared parameters. This sharing of information across 



Table 2 

Posterior means and symmetric 90% credible intervals (CI) for the background heat flows 

qoj (mW/m 2 ) and the mean heat flows vo and us (mW/m 2 ) 





San Rafael Desert 




San Rafael Swell 


Borehole 


Mean 


90% CI 


Borehole 


Mean 


90% CI 


SRD-1 


55.91 


(55.51, 56.32) 


SRS-3 


51.99 


(51.45, 52.54) 


SRD-2 


46.95 


(46.02, 47.90) 


SRS-4 


57.25 


(57.02, 57.47) 


SRD-3 


45.19 


(44.42, 45.92) 


SRS-5 


74.85 


(74.58, 75.11) 


SRD-4 


50.28 


(49.54, 51.01) 


WSR-1 


68.13 


(67.97, 68.29) 


SRD-7 


45.16 


(44.69, 45.63) 








vn 


55.95 


(31.94, 80.39) 


vs 


58.05 


(32.93, 83.08) 
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Fig. 7. Estimated posterior densities of the means, vo and vs (left), and the standard 
deviations, td and ts (right), of the background heat flow for both regions. The (marginal) 
prior distributions are the same for both regions and are presented with dotted lines. 

groups of data is often called borrowing strength and can often lead to bet- 
ter parameter estimates. However, the question here becomes how much 
strength, if any, is borrowed between the nine boreholes. To assess this as- 
pect of the model, we fit single-site models described in Section 3.2 to each 
of the 9 boreholes. 

By performing separate single-site models, we of course do not model 
parameters as spatially dependent. Operationally, some parameters treated 
as random in the combined analysis are assigned fixed values. For example, 
we assign values to \x and T in the prior for the GST histories [see (15)], 
as compared to the additional stage involving fi D , fi s , 7^ and 7! [see (25) 
and (28)] in the combined analysis. To make the results comparable, we set 
the prior mean and covariance matrix of the GST histories equal to their 
marginal prior mean and covariance implied by the multiple-site model. For 
boreholes in the SR Desert (j = 1, . . . , 5), we have the following: 

(29) E(T hj ) = E(E(T hj \iA D )) = E(i* D ) =fx = 0, 

Cov(T^) = Cov(E(T hj \n D )) + E(Cov(T hj \v D )) 

(30) = Cov(/x D ) + E(rtll K ) = (a 2 D + al)I K + E(j 2 D )I K 
= (0.2 + 0.1 + 0.8)I K = 1.1I K . 

For boreholes in the SR Swell (j = 6, . . . , 9), we also have E(Thj) = E(fi s ) = 
and Cov(T/y) = 1.11k- Hence, we set the following prior for T^ in the 
single-site models (same for every borehole): 

(31) T h ~N{n,T)=N K {Q,l.lI K ). 
Similarly, we set the following prior for each q$: 



(32) 



7V(0.06, 0.02 2 + 0.01 2 + 0.1 2 = 0.0105). 
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Finally, we used the same priors for measurement and model error variances 
[see (13) and (14)] as for the multiple-site model. 

We fitted the single-site models via Gibbs samplers, obtaining 10,000 
MCMC samples from the posterior distribution for each borehole and deleted 
2000 iterations for burn- in. 

The estimated posterior means and credible sets for the GST histories 
from both the multiple-site and the 9 single-site models are shown in Fig- 
ure 8. Focusing on the five boreholes in the SR Desert, we see two ma- 
jor differences. First, the GST posterior means (white squares in Figure 8) 
are slightly more dampened for the multiple-site model than the single-site 
models. In other words, we see shrinkage in the posterior means when we 
combine the boreholes. Second, as indicated by narrower 50% and 90% pos- 
terior credible intervals, the posterior uncertainty is substantially less for 
the multiple-site model than the single-site models. We conclude that by 
combining the boreholes in the SR Desert, the GST history parameters are 
borrowing strength across boreholes. 

In the SR Swell the posterior results are quite similar for the multiple-site 
model and the single-site models. In particular, the posterior uncertainties 
are similar in that the 50% and 90% posterior credible intervals are only 
slightly wider for the single-site models than the multiple-site model. We 
believe that the reason for these different results for the two regions is the 
following: The SR Desert GST histories are comparatively similar, so the 
analysis amplifies combining and borrowing strength. The SR Swell results 
seem more diverse, suggesting borrowing strength should be comparatively 
weak. 

Posterior density estimates (not shown here) for the background heat 
flows qoj and the variances Oy and a 2 were almost identical for the two 
analyses. 

4.3. Sensitivity analyses. In our analyses we used prespecified fixed val- 
ues of the temperature intercepts Tqj and hyperparameters (i.e., fixed pa- 
rameters of prior distributions). The selection of hyperparameters was dis- 
cussed in Section 3.4. We next assess the sensitivity of results to some of 
these specifications. 

4.3.1. Temperature intercept Tqj. The temperature intercepts Tqj used 
here were least square estimates of the intercepts in simple linear regressions. 
That is, separately for each borehole, temperature data were regressed on 
the thermal resistance vector Rj. In these steps, the regressions were based 
only on data at the deep parts of the borehole, specifically below 150 meters 
for boreholes in the SR Desert and below 200 meters for boreholes in the 
SR Swell. We used the usual least squares standard errors to guide our 
sensitivity analysis on Tqj. Specifically, we fitted the multiple-site model in 
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Fig. 8. Comparison of posterior distributions of GST histories from the multiple-site 
model (black) and single-site models (grey). The white squares show the estimated posterior 
means and the vertical bars show the symmetric 50% (thick) and 90% (thin) posterior 
credible intervals. 



two additional cases: (1) all Toj's were set to three standard errors below the 
least squares estimates and (2) all Trjj's set to three standard errors above 
the least squares estimates. 

Overall the results were not highly sensitive to these changes in the tem- 
perature intercepts. Estimated densities of the heat flow parameters q^j (not 
shown here) indicated that posterior for the qoj responded to changes in the 
Tqj's as we expect a slope to change when the intercept is changed. When 
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Fig. 9. Comparison of posterior distributions of GST histories for 4 boreholes using 
different values for the temperature intercept Try . The white squares show the estimated 
posterior means and the vertical bars show the symmetric 50% (thick) and 90% (thin) 
posterior credible intervals. At each time point the middle bar shows the original results 
(same as in Figure 2) and the left and right bars show the results for lower and higher 
values of Tiy , respectively. 

the Tqj were lowered, the qoj were higher and vice versa. However, posterior 
means and standard deviations of the heat flows, ud, v$, td, ts, and ayj 
and <7j, did not vary much as we changed Tqj's. 

The estimated posterior means and credible sets for the GST histories 
for the three cases of Tqj are shown in Figure 9. We show results for two 
boreholes from each region, but the effects are similar for the other boreholes. 
The main effect of changing the Trjj's is that the posterior means of the 
GST histories are slightly shifted. The temperature anomalies are higher 
when Tqj is lower and vice versa, so it seems that they are compensating 
the changing Tqj. Interestingly, the amount of shifting increases as we go 
further back in time. On the other hand, the changes in posterior means 
are not large when compared to the posterior uncertainty of the results. We 
conclude that sensitivities to the specifications of Tqj are overshadowed by 
the posterior uncertainty of the results. 



4.3.2. Hyperparameters i]o, rjs> Vo- Recalling (24), the regional prior 
means, vd and vs, of the heat flows have prior standard deviations (r?^ 



^2)0.5 

?/i)) a5 



and (Vs + Vq) '' 



respectively, and prior correlation ^/((^o + 7 ?i>)( r ?o + 
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Table 3 

The original and four new settings of the hyperparameters rf^,, rj s and n% along with the 

implied prior standard deviations of vd and vs and the prior correlations between vd 

and vs 





vl 


and 


ri% 


vl 


Prior sd. [mW/m 2 ] 


Prior cor 


Original setting 




10 2 




20 2 


22.36 


0.80 


Setting 1 




20 2 




20 2 


28.28 


0.50 


Setting 2 




30 2 




20 2 


36.06 


0.31 


Setting 3 




100 2 




O 2 


100.00 


0.00 


Setting 4 




100 2 




67- 


120.37 


0.31 



We report on results for four additional settings of 770, r/D and rjs leading 
to the prior standard deviations and correlations given in Table 3. Note 
that in light of the range of heat flow estimates found in the literature (see 
Section 3.4), a prior standard deviation of 120 mW/m 2 is very large. Also, 
the prior correlation we selected for the original setting is very high (0.8). 
Note also that sensitivity settings 2 and 4 have the same correlation but 
different standard deviations. 

First, posterior means and credible intervals of the GST histories (T^,-; 
not shown) were almost identical across all five cases. The same was true 
for the means and standard deviations of the histories ((/>£>, Hd-> Id and 75; 
not shown) Also, the posterior density estimates (not shown) for heat flow 
parameters qj were almost identical in all cases. 

The posterior distributions of vp) and v$ displayed strong sensitivities. 
Figure 10 (right) shows the posterior means and credible intervals for vd 
and vs for all 5 settings. Not surprisingly, as the prior uncertainty increases, 
the posterior uncertainty increases. We note differences in the posterior 
means that correspond to the different prior correlations. In the original 
setting the posterior means of vp> and vs were very close and seemed not 
to respond to the regional information indicated in the posteriors of the qoj 
[the vertical line segments in Figure 10 (right) are the posterior means of 
the qoj] also see Table 2]. Note that as the prior correlation decreases, the 
posterior means of vp) and v$ separate and approach the averages of the heat 
flow posterior means in their corresponding regions. However, due to the rel- 
atively large posterior uncertainties of vp> and v$ in all cases, the changes 
we see in the posterior means are all within the 50% credible interval of the 
original model. 

While these results are not surprising, it is instructive to see the workings 
of Bayesian updating of borrowing-strength priors. Finally, though the sam- 
ple sizes at each borehole are relatively large, the operative "sample sizes" 
for treating vp> and vs are 5 and 4 sites, respectively. 
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Fig. 10. Left: Comparison of posterior distributions of the mean GST histories, fj, D 
and fj, s , using different values of a D , ag and o~q (see Table 4). The white squares show the 
estimated posterior means and the vertical bars show the symmetric 50% (thick) and 90% 
(thin) posterior credible intervals. Right: Comparison of mean heat flow vo and vs using 
different values of 7)%, rfg and rfe (see Table 3). The white squares show the estimated 
posterior means and the vertical bars show the symmetric 50% (thick) and 90% (thin) 
posterior credible intervals. The line segments at the bottom show the posterior means of 
the heat flows qoj . 



4.3.3. Hyperparameters or>, as, 0q. The hyperparameters <td, &s and ctq 



,2\0.5 

>o) > 



of 



determine the prior standard deviations, (p D + ctq) and (a s + a, 
the elements of the prior mean GST history vectors Hd and /15 . They also 
imply that the prior correlations con (fiDk^Sk), k = 1, . . . ,K, are equal to 
Oq/({oq + cr |?)( cr o + <T I))°' 5 - We considered three additional settings for these 
parameters as shown in Table 4. 

Overall, the results were not sensitive to the different settings of <td, 175 
and <7q. The mean GST histories in each region, fj, D and n s , were only 



Table 4 

The original setting of the hyperparameters a 2 D , a% and oq and the three additional 

settings along with the implied prior standard deviations of the elements of the vectors 

fj, D and fj, s and the prior correlations between fiok and fisk for k — 1, . . . , K 



<T 2 



Prior sd. (°C) 



Prior cor 



Original setting 


0.2 


0.1 


0.55 


0.33 


Setting 1 


0.1 


0.1 


0.45 


0.50 


Setting 2 


0.2 


0.0 


0.45 


0.00 


Setting 3 


0.3 


0.15 


0.67 


0.33 
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slightly affected by the changes in the prior [see Figure 10 (left)]. The poste- 
rior credible intervals are slightly wider when the prior standard deviations 
are higher. The posterior means are almost identical for the four settings 
over the first three or four centuries. In the last century there is a slight dif- 
ference between settings 1 and 2, particularly in the SR Swell. These settings 
have the same prior standard deviation but different prior correlations (0.5 
and 0, resp.). It seems that higher prior correlation leads to posteriors that 
favor similarity of the elements of fi D and /jg , but only for the more recent 
time points. Again, we note that these differences are very small compared 
to the posterior uncertainty in the results. 

The results for the GST histories (T/y, not shown here) were very similar 
for all four settings. The only differences were that the posterior uncertain- 
ties increased slightly when the prior standard deviations of [i D and fi s 
increased. The different prior correlations did not seem to matter, settings 1 
and 2 gave almost identical results. Density estimates (not shown) for other 
parameters were virtually identical for the different prior settings. 

4.3.4. Correlated measurement and model errors. Though we expect sub- 
stantial correlation among the elements of each Y,- and among the elements 
of each T r j , we assumed conditionally independent measurement errors and 
model errors in (16) and (17), respectively. The modeling notions are that 
the lion's shares of these correlations are due to the structure of the true tem- 
peratures in the case of the Yj and the structure in temperatures captured 
by the heat equation model in the case of the T r j [also see Section 3.2(iii)]. 
Neither notion is unassailable: while the errors in (16) are due to the mea- 
surement process, they also respond to approximation errors associated with 
the use of the reduced temperature definition. Similarly, the errors in (17) 
are attributable to model approximations associated with the specifics of our 
use of the heat equation. However, we have virtually no prior information 
regarding the structure of the unmodeled physical processes; if we did know 
more, that knowledge could be used to improve the physical models. 

Ignorance is not a valid defense for independence assumptions. Rather, 
those assumptions lead to obvious simplifications in the analysis and avoid 
difficulties associated with potential nonidentifiability issues. However, some 
sensitivity checks are desirable. In our setting the major concern is that in- 
correct assumptions of independence may lead to underestimation of uncer- 
tainties in the final results. 

To assess independence assumptions, we examined model "residuals." 
Some indication of structure may be developed by inspecting estimated er- 
rors defined by 

(33) e, = Yj - [E(T rj ) + T 0j 1 N . + E(q oj )Rj], j = 1, . . . , 9, 
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using estimated posterior expectations as indicated. A more appropriate 
approach is to compute 

(34) e™ = Y, - [T™ + T oj l Nj + q^R,], j = 1, . . . , 9, 

where the superscripts m index MCMC iterations, though this option leads 
to an ensemble of residuals. 

For each of the nine boreholes, we fitted time-series style ARMA models 
to the 6j. Though we noted some differences, AR(1) and ARMA(1,1) models 
provided reasonable fits. Since the forms of the covariance matrices of the 
AR(1) and ARMA(1,1) are quite similar, we focused on AR(1) models. We 
also inspected realizations of residuals as defined in (34). These residuals 
were similar to those based on (33), though they exhibited less structure, 
suggesting that basing our sensitivity checks on (33) is conservative. We 
replaced the covariances in (17) by the covariance matrices cr?C(0), where 
C(4>) is a correlation matrix with ones on the diagonal and off diagonal ele- 
ments (f) k for every pair of depths k units apart (here one unit corresponds to 
5 meters). We treated <j> as a known quantity and reran the MCMC analyses 
for two choices of <p: 0.65 and 0.85. The final posterior results were not very 
different from those using the independence assumption. In Figure 11 we 
display comparisons of posterior distributions of GST histories for two bore- 
holes for each of the Desert and Swell regions as well as the corresponding 
region mean processes. The boreholes were selected to indicate the range of 
differences, that is, one borehole showing the least differences and another 
one showing the most differences. We did observe nonnegligible increases in 
the spreads of the posteriors for the heat flow parameters qoj, though there 
were virtually no changes in their means. 

We repeated this process by replacing the model error covariances in (16) 
by matrices a\. C(^) for the same two values of <fi. The same behaviors 
were observed; indeed, the differences were smaller that those above. Hence, 
though some structure in the errors are unmodeled in this article, the effects 
of this appear to be very minor in the posterior inferences regarding GST 
histories. 

5. Discussion. We developed Bayesian hierarchical models featuring two 
key aspects: (1) the use of a physics based model having uncertain parame- 
ters and subject to model error, and (2) the illustration of borrowing strength 
based on priors on selected parameters to combine data sets. In the first 
development we relied on a common framework to define reduced temper- 
atures T r to justify use of the heat equation as the main physical model. 
However, we added the treatment of background surface heat flows go as un- 
known parameters and the inclusion of model errors implying a random heat 
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Fig. 11. Comparison of posterior distributions of GST histories for 4 boreholes and the 
mean histories fi D and fi s using model error covariance matrices cr?C(^>) with different 
values of cj>. The white squares show the estimated posterior means and the vertical bars 
show the symmetric 50% (thick) and 90% (thin) posterior credible intervals. At each time 
point the first (darkest) bar shows the original results (<j> — 0; same as in Figure 2) and 
the next two bars show the results for <f> = 0.65 and <f> — 0.85. 



equation model. The hierarchical modeling approach also allows us to sepa- 
rate explicitly measurement errors from model errors. These steps support 
the claim that our models account for important uncertainties. 

In our prior distributions we modeled the surface history processes and 
background heat flows as arising from region-specific (SR Desert or SR Swell) 
models, which in turn have parameters generated from a San Rafael basin- 
wide prior model. These formulations led to posterior results that display 
very notable and appropriate behaviors. As discussed in Section 4.2, the in- 
ferences for model parameters and histories in the SR Desert region indicated 
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borrowing strength in concert with the similarities of behaviors at individual 
boreholes. By comparison, we noted very little borrowing of strength in the 
SR Swell region where the individual results were comparatively dissimilar. 
We find this result quite satisfactory and illustrative, but should note that 
the priors used were extremely simple. They were chosen because we had 
little prior information and too little data (i.e., four and five boreholes in the 
two regions) on which to update more intense priors. When feasible, we rec- 
ommend consideration of more intense spatial process priors for parameters. 
See Cressie (1993) and Banerjee, Carlin and Gelfand (2004) for discussion 
of sophisticated spatial models. 

As is appropriate in most Bayesian analyses, we devoted substantial at- 
tention to sensitivity analyses. The results are discussed in Section 4 and not 
reviewed here. However, note that we did not present analyses regarding the 
specification of the thermal conductivities (k) used in defining the vectors 
of thermal resistance vectors R. Very cursory inspections suggest to us that 
the approach may be very sensitive to these quantities. We will pursue this 
aspect in an alternative modeling approach to be reported on elsewhere. 

We note that the results lead to the suggestion that borehole data are 
useful in inferring surface temperatures for times from the recent past to 
about 200 years in the past. For times deeper in the past, the borehole 
data appear to be comparatively less informative. We base this suggestion 
on the behavior of the posterior distributions of the site-wise histories as 
well as the SR Desert and Swell mean histories. These distributions seem 
to asymptote to region-specific distributions. For all times before 1800 and 
all five boreholes in the SR Desert, the posterior standard deviations of 
historical temperatures vary between 0.61 and 0.68 (only 5 of the 35 values 
are less than 0.65); these values are 2 or 3 times larger than the standard 
deviations for temperatures in the Desert in 1980. In the SR Swell, all 28 
of the corresponding standard deviations are between 0.87 and 0.95, and 
are roughly 4 times the standard deviations for temperatures in the Swell in 
1980. Regarding this issue, North et al. [(2006), page 80] write the following: 

The time resolution and length of borehole-based surface temperature recon- 
structions are severely limited by the physics of the heat transfer process. . . . 
A surface temperature signal is irrecoverably smeared as it is transferred to 
depth. The time resolution of the reconstruction thus decreases backward in 
time. For rock and permafrost boreholes, this resolution is a few decades at 
the start of the 20th century and a few centuries at 1500. 

For related discussion see Beltrami and Mareschal (1995) and Hopcroft, 
Gallagher and Pain (2007). Our addition to these claims is that the posterior 
standard deviations based on models that include model error and other 
uncertainties are roughly constant by necessity for times beyond 200 years 
in the past. Of course, this is based on limited data from a limited region 
and need not apply in greater generality. 
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To characterize our results in regard to climate change, we provide in- 
ferences on the changes in surface temperatures over the four periods 1600, 
1700, 1800 and 1900 to the latest year in our data sets (1980) for each of the 
nine boreholes and for the SR Desert and SR Swell means (see Figure 12). 
We note that the point estimates of these changes are typically positive 
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Fig. 12. Posterior means and credible sets of surface temperature changes over the four 
periods 1600, 1700, 1800 and 1900 to the latest year in our data sets (1980) for each of 
the nine boreholes and for the SR Desert and SR Swell means. The vertical bars show 
the symmetric 50% (thicker) and 90% (thinner) posterior credible intervals and the grey 
horizontal lines show the posterior means. 
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(i.e., increased temperature) for the individual boreholes and are all posi- 
tive for the basin-wide means. However, all 90% credible intervals cover 0°C, 
though many of the 50% credible intervals lie above 0°C. Based on the com- 
mon trend in these results, we believe that the suggestion of warming is 
supported. While the strength of this support is not strong, we note that 
our sample sizes are very small. Further, since our data ends in 1980, we 
cannot find the more recent warming reflected in other data. Finally, we 
caution that traditional quantification associated with so-called statistical 
significance (i.e., 90% or 95% intervals not covering 0) are of little relevance 
in regard to decision making in the context of climate change. 

APPENDIX A: LIST OF PARAMETERS 

Unknown parameters: 

Borehole specific parameters (j = 1, . . . , 9): 

Thj'- Temperature histories, vectors of dimension K. 
T rj -: True reduced temperatures, vectors of dimension Nj. 

qoj: Heat flows, scalars. 
o~y~'- Measurement error variances. 

a?: Model error variances. 

Region specific parameters: 

/x D , fi s : Mean temperature histories for SR Desert (D) and SR Swell (S), 

vectors of dimension K. 
7z>> Ts : Variances of temperature histories for SR Desert (D) and SR Swell 

(S). 
up, v$: Mean heat flow for SR Desert (D) and Swell (S), scalars. 
r D' T s' Variances of heat flow for SR Desert (D) and Swell (S). 

Hyperparameters-constants: 

Prior mean of /x D and n s . 

Define the prior covariance structure of fi D and /i s . 

Prior mean of vd and v$. 

Define the prior covariance structure of vjj and v$. 

Define the Inverse Gamma prior for a Y 7 -- 

Define the Inverse Gamma prior for ct? 

Define the Inverse Gamma prior for 7^ and 7I. 

Define the Inverse Gamma prior for r|> and r|. 

APPENDIX B: GIBBS SAMPLER 

A Gibbs sampler is a method that obtains approximate samples from 
the posterior distribution. It avoids the big dimensionality of the model 
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by simulating only parts of the parameters at a time, using the so-called 
full conditional distributions. The notation [X\-] reads "the full conditional 
distribution of X given all other parameters and the data." Also, ||x|| = x'x, 
for a vector x. 

GST history vectors T^j and reduced temperature vectors T r j. To 
make the Gibbs sampler more efficient, we sample the joint full conditional 
distribution Thj,T r j\- for each borehole j. Note that the joint (prior) dis- 
tribution of Tf L j and T r j for boreholes in the SR desert (J = 1, . . . , 5) is 

<» ( T <)~<(t:)lh it) 

where 

(36) E j = a]l Nj + 1 2 D A j A' j . 

The joint full conditional distribution of T r j and Tfij given Yj and all other 
parameters, which we denote by G, can be written as follows: 

\t t iv ial — L_i I zll ' _J \~Ii I ' it it ai 
[i- rj ,l hj \Yj,(d>\- [L hj \L rj ,V\ 

(37) lY ' |OJ 

— [Trj | Yj , 0] [T h j | T r j , 0] . 

Therefore, to sample T/y,T r j|-, we first sample the marginal distribution 
[TVj|Yj,0] and then the conditional distribution [T/y|T r j,0]. First note 
that [T rj |Y,-,0] is proportional to [Yj \ T r j , 9] [T rj \ 0] where [T rj |9] is the 
marginal distribution of the joint prior (35). Therefore, for boreholes in the 
SR Desert (j = 1, . . . , 5), we have [T rj \Yj,Q] = N(Dd, D), where 

/ 1 ~ ^- x 

(38) D = -f-I N . + EJ 

and 

(39) d = (Yj - (T 0j l Nj + qjRj))/^ + EtU,-/^. 
Second, note that [T/y|T r j,0] is the conditional prior distribution 

(40) N(» D + -ftAfiftTrj - A^ D ) n ll K - -ybAfi^Aj). 

For boreholes in the SR Swell (j = 6, . . . , 9), we replace \l d in (39) and (40) 
with fi s and 7^ in (40) and (36) with 7^. 

Measurement and model error variances, er^ ■ and a?. 

(41) a Yj \- ~Ig(J^ + a Y ,b Y + \\\Yj - (T rj + T oj l Nj + <&R,-)IM > 

(42) otl-nlG^ + aJ+^llTrj-AjThjU, j = l,...,9. 
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Mean and variances of the GST history, fj, D , fj, s , 7^ and ~/p. We sample 
the joint full conditional distribution of fi D and fi s . We have fj, D ,fi s \- ~ 
N(Dd,D) with 

(43) 
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where ^ = (^ + ^)(^ + ^)-^ and W = K + a 2 )(4 + a 2 )-a 4 . 
The variances of the GST histories are sampled separately for each region, 

(45) 
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and 
(46) 



/ 1 9 
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Heat flow parameters goj? ^Di ^s? t £)' t S' ^or boreholes in the SR 
Desert (J = 1, . . . , 5), we have 



(47) gj \-~N( D jK 3 



■ r.J 



T 0j Inj ) + o-yj ^d Td^Yj 



r|,R^ + 4. 



'^R-iRi + 4^ 



For boreholes in the SR Swell (j = 6, . . . ,9), we replace fj, D and t^ in (47) 
by fi s and r|. 

The mean heat flow parameters vd and us are sampled from a joint 
distribution I'd, 1/5 ]■ ~ N(Dd,D) with 



/ 4 r? 2 + r,l 



(48) 
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and 



(49) 
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T D j=1 
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where w = (£ + E±3i)(^ + 1±2l) _ % and « = (r? 2 + r/ 2 ) 2 

Finally, the variances of the heat flow are sampled separately for each 
region, 



(50) 

and 

(51) 
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